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ABSTRACT 

A  key  aspect  of  estimating  transmissivity  distributions  from  pumping  tests  is  the  spatial  resolution  of  the 
estimation.  Some  assumption  of  spatial  structure  is  needed  to  make  an  inversion  possible.  Both  the  nature 
of  groundwater  flow  and  assumed  structure  cause  a  smeared  estimate.  This  research  is  focused  on  the 
how  information  content  of  the  data  smears  the  estimate.  The  inversion  procedure  near  the  solution  is 
approximated  as  iterative  linear  updates.  Linear  resolution  analysis  of  the  distributed  sensitivity  to 
transmissivity  indicates  the  limit  of  detail  possible  from  the  inversion. 

To  enable  this  research,  adjoint-sensitivity  calculations  were  added  to  MODFLOW-2000.  The 
adjoint-sensitivity  calculations  significantly  speed  the  calculation  of  the  sensitivity  matrix. 

To  highlight  the  concepts,  we  use  a  homogeneous  aquifer  with  known  storativity.  The  analysis  indicates 
resolution  of  transmissivity  is  finer  near  the  pumping  and  monitoring  wells  than  in  the  regions  between 
wells.  Transient  data  are  highly  redundant  with  respect  to  transmissivity  information.  Very  early  time  data 
do  provide  different  spatial  information  than  later  time  data  do. 

INTRODUCTION 

An  important  question  in  any  endeavor  is  “How  well  can  we  do  it?”  This  is  often  a  very  complicated 
question.  Bounding  the  answer  rather  than  answering  it  directly  often  turns  a  difficult  question  into  one  we 
can  answer.  Our  endeavor  is  to  estimate  subsurface  properties.  Resolution  analysis  provides  insight  into 
how  well  we  can  estimate  the  properties.  Resolution  analysis  is  a  common  technique  in  the  geophysical 
community  (Menke,  1989;  Tarantola,  1987).  With  the  exception  of  Vasco  et.  at.  (1997),  resolution  analysis 
has  not  received  attention  in  the  groundwater  literature. 

Resolution  analysis  uses  a  linear  approximation  of  the  dependence  of  measurements  on  estimated 
parameters  to  determine  how  the  inversion  estimates  relate  to  the  actual  subsurface  properties.  This  is 
truly  a  bounding  analysis  because  we  assume  that  our  groundwater  model  is  accurate  and  that  the 
estimated  parameters  are  close  enough  to  correct  that  non-linearity  can  be  ignored.  We  also  partially 
neglect  the  influence  of  errors  in  the  measurements  used  in  the  inversion. 

In  this  paper,  we  consider  the  estimation  of  the  distribution  of  transmissivity  in  a  hypothetical  aquifer  from 
drawdown  measurements  in  observation  wells  about  a  pumping  well.  The  emphasis  is  on  the  resolution 
analysis  not  on  a  specific  inversion  procedure.  As  this  is  a  conference  on  MODFLOW,  the  addition  of  the 
adjoint  state  method  (Sun,  1994;  Townley  and  Wilson,  1985)  of  calculating  sensitives  to  MODFLOW-2000 
(Harbaugh  et  al.,  2000;  Hill  et  al.,  2000)  and  its  use  in  this  investigation  is  highlighted. 

SENSITIVITY  CALCULATIONS 

The  simulated  hypothetical  aquifer  has  a  central  region  12.9  m  by  7.6  m  of  constant  grid  spacing  of  0.25 
m.  Twelve  wells  are  positioned  in  roughly  concentric  rings  about  the  central  well  as  shown  in  Figure  1 .  The 
central  region  is  surrounded  by  a  grid  that  expands  with  a  growth  ratio  of  1 .5  to  1 50  m.  The  grid  extends 
beyond  the  central  region  far  enough  that  no  appreciable  drawdown  occurs  at  the  boundary  during  the 
simulation  period.  A  constant  storativity  of  0.3  is  used.  A  uniform  transmissivity  of  0.0032  m2  /s  is  used.  A 
constant  pumping  rate  of  5  I/s  is  maintained  for  10,000  s  (170  min).  No  appreciable  drawdown  occurs  at 
the  boundary  during  the  simulation  period.  A  total  of  122  observations  were  made  from  drawdown  at  the 
wells.  Four  drawdown  observations  per  log  base  ten  time  interval  were  recorded.  Early  observations, 
before  a  well  responds,  were  not  included. 
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The  model  has  a  total  of  9409  cells.  Estimating  a  distribution  of 
transmissivity  refers  to  allowing  each  cell  of  the  MODFLOW  model 
to  be  assigned  a  different  value.  As  a  consequence,  a  matrix  of 
parameter  sensitivities,  the  sensitivities  of  the  drawdown 
observations  to  transmissivity,  must  be  calculated  for  each  cell  in  the 
model.  Calculating  the  sensitivity  matrix  using  the  sensitivity 
equation  process  of  MODFLOW-2000  requires  a  separate  calculation 
for  each  cell.  This  calculation  required  19  h  on  a  1  GHz  PC.  These 
calculations  agree  well  with  the  analytic  solutions  presented  by  Oliver 
(1993).  An  adjoint  states  based  calculation  of  drawdown  sensitivity  to 
hydraulic  conductivity  was  added  to  MODFLOW-2000  following  the 
procedure  described  in  Carrera  and  Medina  (1994). 

The  procedure  of  Carrera  and  Medina  requires  one  adjoint  states 
calculation  per  observation  location.  These  calculations  took  2.3  min¬ 
utes.  The  calculations  were  not  sufficiently  accurate  for  our  purpose. 
They  disagree  by  10%  to  20%  with  the  sensitivity  equation  based 
calculations.  Errors  were  larger  for  cells  with  small  sensitivity  but  these  errors  are  not  relevant.  The 
inaccuracy  may  represent  a  coding  error  or  a  fundamental  limitation.  Improvements  to  the  calculations  are 
being  investigated.  While  not  sufficient  for  the  resolution  analysis,  20%  accuracy  may  be  acceptable  for  the 
early  stages  of  an  inversion  where  non-linearity  makes  any  calculation  inaccurate. 

The  adjoint  sensitivity  calculations  that  were  used  were  based  on  a  separate  adjoint  states  calculation  for 
each  of  the  122  observations.  These  calculations  took  94  minutes  to  complete,  12  times  faster  than  the 
sensitivity  equation  calculation.  Errors  for  significant  sensitives  were  less  than  0.1%  for  observations  that 
coincided  with  the  end  of  a  time  step  and  between  0.1%  and  2%  for  interpolated  observations. 
Improvements  to  speed  are  anticipated.  We  have  noted  that  if  only  the  first  third  of  the  adjoint  states 
calculation  were  performed  an  addition  factor  of  three  speed-up  could  be  obtained  without  appreciable 
increase  in  error. 
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Figure  1  Well  Field 
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RESOLUTION  MATRIX 

The  resolution  matrix  is  based  on  the  sensitivity  calculation  of  the  last  step  of  a  non-linear  inversion  proce¬ 
dure  (Tarantola,  1987).  At  the  last  step  the  inversion  can  be  linearized  to 


Gm  =  d  (1) 

where  G  is  the  sensitivity  matrix,  m  is  an  unknown  error  in  the  transmissivity  parameters,  and  d  is  the 
difference  between  the  observation  with  the  current  parameter  estimates  and  the  measured  observations. 
If  G  could  be  inverted  ,we  could  update  the  transmissivity  estimate  by  Am  =  G_1d,  but  G  is  a  138  by  9409 
matrix  which  does  not  have  an  inverse.  Instead,  we  write  the  update  as 


Am  =  G  sd 


(2) 


G  9  is  referred  to  as  the  generalized  inverse.  The  update,  Am,  can  be  related  to  the  error,  m,  by 
combining  1  and  2. 


Am  =  G  sGm  (3) 

R  =  G~SG  is  the  resolution  matrix  (Menke,  1 989;  Tarantola,  1 987).  The  resolution  matrix  describes  how  the 
transmissivity  update  is  related  to  the  transmissivity  error.  Because  G~s  is  not  the  inverse  of  G,  the  update 
is  not  equal  to  the  error.  It  can  be  shown  that,  if  the  transmissivity  estimate  is  nearly  correct  then  non-linear 
effects  are  not  important  and  equation  3  holds  for  the  transmissivity  estimate  itself.  That  is, 

Test  =  RTtrue  (4) 
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where  Test  is  the  transmissivity  estimate  of  the  inversion  and  Ttrueis  the  real  transmissivity.  Equation  4 
neglects  errors  in  the  model,  measurements  errors,  and  other  plagues  of  inversion.  These  caveats  aside, 
the  resolution  matrix  provides  us  with  a  tool  that  we  can  use  to  set  our  expectations,  design  experiments, 
and  compare  different  inversion  strategies. 

ANALYSIS  USING  TRUNCATED  SINGULAR  VALUE  DECOMPOSITION 

In  equation  2,  G~g  incorporates  the  constraints  used  in  the  inversion.  These 
constraints  may  include  prior  information,  regularization  equations  to  impose  a 
structure  on  the  solution,  or  a  prior  covariance  structure  arising  from  a  Bayesian 
formulation  of  the  inversion  problem.  The  constraints  must  be  included  in  the 
resolution  analysis  of  a  particular  inversion  process.  We  simplify  the  analysis 
by  using  truncated  singular  value  decomposition  as  the  way  of  calculating  the 
generalized  inverse.  Menke  (1 989)  refers  to  this  G~s  as  the  natural  generalized 
inverse.  It  is  natural  because  it  constrains  the  estimated  distribution  to  only 
have  variations  that  are  distinguishable  from  the  measurements.  Any  form  of 
regularization  introduces  a  bias  into  the  estimation.  The  bias  of  the  truncated 
singular  decomposition  approach  is  based  on  the  information  content  in  the 
data.  If  a  variation  in  transmissivity  cannot  be  clearly  supported  by  the  data,  it 
is  left  out  of  the  solution. 

The  singular  value  decomposition  of  G  results  in 

G  =  UAVT  (5) 

where  A  is  a  diagonal  matrix  of  eigenvalues,  A;,  and  U  and  V  are  matrices  containing  the  associated 
eigenvectors.  We  refer  to  U  as  the  data  eigenvectors  and  V  as  the  transmissivity  eigenvectors. 

Because  of  the  orthogonality  and  normality  properties  of  U  and  V 

G-g  =  VA-1UT  (6) 

where  A-1  is  a  diagonal  matrix  with  diagonal  elements  j-.  Very  small  eigenvalues  indicate  a  redundancy 
in  the  information  about  the  transmissivity  contained  in  the  observations.  Two  measurements  with  exactly 
the  same  information  content  would  produce  an  eigenvalue  with  a  value  of  zero. 

Figure  2  plots  the  normalized  eigenvalues  from  singular  decomposition  of  the  sensitivity  matrix  using  a 
semi-log  scale.  The  eigenvectors  are  normalized  with  respect  to  the  largest  eigenvalue.  The  figure  shows 
a  rapid  decline  in  the  eigenvalues  implying  a  large  redundancy  in  information.  Two  sets  of  eigenvalues  are 
plotted.  One  using  all  the  138  observations  and  one  using  only  the  13  observations  acquired  at  10,000  s, 
the  last  time  step  of  the  simulation.  The  eigenvalues  of  these  13  observations  nearly  coincide  with  the 
thirteen  largest  eigenvalues  of  the  full  data  set.  The  overlap  of  the  two  sets  suggests  that  observations 
from  different  locations  provide  much  more  new  information  than  can  be  obtained  by  adding  a  large 
number  of  new  observations  at  current  locations,  a  conclusion  that  is  supported  by  the  resolution  analysis 
presented  below. 

A  consequence  of  small  eigenvalues  is  that  there  are  weighted  combinations  of  transmissivity,  Am,  that 
have  a  very  small  influence  on  the  calculated  drawdown  observations.  These  combinations  are  defined  by 
aVi,  where  a  is  a  constant.  There  is  a  corresponding  weighted  combination  of  drawdown  measurements, 
d],  defined  by  /?U;  that  are  influenced  by  the  V;.  If  d;  is  imperfect,  as  of  course  it  will  be  using  measured 
data,  then  the  errors  in  d;  will  be  magnified  by  in  the  resulting  prediction  of  Am;.  To  avoid  this  problem, 
we  truncate  the  singular  value  decomposition  at  some  minimum  A].  This  is  the  method  whereby  only 
variations  in  transmissivity  that  can  be  supported  by  the  data  are  retained.  The  appropriate  truncation  level 
is  related  to  the  uncertainty  in  the  observations.  Vogel  (2003)  provides  an  algorithm  for  selecting  an 
appropriate  truncation  level. 
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Figure  2  Normalized 
eigenvalues  obtained 
from  singular 
decomposition  of  the 
sensitivity  matrix. 
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RESOLUTION  OF  TRANSMISSIVITY  ESTIMATES 

The  resolution  matrix  is  large;  number  of  parameters  by  number  of  parameters.  The  spread  function 
provides  one  way  of  reducing  the  resolution  matrix  to  a  comprehensible  level.  Menke(  1989,  pg  68)  defines 
the  spread  function  for  the  resolution  matrix  as 


spread^^^Rij-lij)2  (7) 

i=i  j=i 

where  I  is  the  identity  matrix.  We  define  the  spread  of  each  cell  as 

n 

Si  =  1  -  ^  (Rij  -  lij)2  (8) 

j=i 

A  value  of  1  for  a  cell  implies  that  the  cell  can  be  estimated  with  good  accuracy.  It  also  implies  that  the 
estimate  for  the  cell  is  not  influenced  by  the  transmissivity  of  other  cells.  A  small  value  implies  that  the 
estimates  of  transmissivity  for  a  cell  are  strongly  influenced  by  other  cells.  A  value  near  zero  indicates  the 
the  transmissivity  estimate  for  the  cell  will  not  related  to  the  actual  transmissivity  of  the  cell. 


Figure  3  plots  the  spread  functions  for  three  separate  cases;  a)  using  all  the  eigenvectors  of  the  138 
observations,  b)  using  the  largest  29  eigenvectors,  and  c)  using  only  the  13  observations  acquired  at  the 
last  time  step  of  the  simulation.  The  plots  are  color  based  and  each  plot  is  rescaled  so  that  the  largest 
values  are  purple  and  and  small  values  are  red.  The  color  scale  is  drawn  to  the  right  of  each  plot.  Cells 
with  values  below  the  cutoff  on  the  color  scale  are  not  drawn. 
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Figure  3  Spread  functions  of  the  resolution  matrix,  a)  Using  all  the  eigenvectors,  b)  Using 
eigenvectors  of  the  29  largest  eigenvalues  c)  Using  only  10,000  s  observations. 
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Plot  3a  would  be  appropriate  if  the  eigenvalues  were  all  approximately  the  same  size  indicating  that  all  of 
the  eigenvectors  could  be  resolved  by  the  inversion  process.  The  plot  shows  values  near  1  for  the  cells 
immediately  adjacent  to  the  pumping  well  and  relatively  large  values  in  the  vicinity  of  the  well,  indicating 
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that  the  transmissivity  of  a  disturbed  zone  could  be  resolved  at  this  grid  size.  Away  from  the  well  the 
spread  function  diminishes  rapidly  with  only  a  3  m  diameter  circle  showing  a  strong  anticipated  correlation 
between  the  estimate  and  the  actual  local  transmissivity.  The  cells  immediately  adjacent  to  the  observation 
wells  would  also  be  correlated  with  their  estimates.  If  observations  were  acquired  more  often,  then  a  plot 
similar  to  plot  3a  would  indicate  higher  resolution,  an  unrealistic  conclusion. 


Plot  3b  is  more  realistic.  Only  the  29  eigenvectors  associated  with  normalized  eigenvalues  larger  than 
1CT3  were  used  to  create  the  spread  function.  The  plot  indicates  that  nowhere  can  transmissivity  be 
resolved  at  the  25  cm  scale  of  the  grid.  With  a  lower  cutoff  of  0.01 ,  the  large  red  area  does  not  indicate  a 
region  of  even  moderate  resolution.  As  with  the  eigenvalue  comparison  of  the  full  observation  set  with  just 
one  observation  at  each  well,  a  comparison  of  plots  3c  to  3b  indicates  that  the  full  set  of  observations  only 
provides  a  small  improvement  as  we  concluded  earlier  from  Figure  2. 


One  might  conclude  that  transient  data  are  not  of  value  in  estimating  a  transmissivity  distribution.  We  do 
not  believe  that  such  a  conclusion  is  justified.  Investigations  of  the  sensitivity  matrix  indicate  that  sensitivity 
is  different  during  the  early  portion  of  the  drawdown  compared  to  the  late  portion  and  that  the  shape  of  the 
sensitivity  field  changes  rapidly  in  this  region.  This  suggests  that  more  measurements  at  early  times  may 
be  valuable.  The  crux  of  the  matter  lies  in  the  signal  to  noise  ratio  of  the  measurements.  A  larger  signal  to 
noise  ratio  correlates  with  a  lower  eigenvalue  cutoff,  more  retained  eigenvectors  and  more  finely 
determined  transmissivity.  Unfortunately,  the  most  important  data  is  related  to  small  values  of  drawdown. 
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“""ive  sensitivity." 
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Weii^/n  presented  resolution  analysis  as  a  tool  for  inv<  istigaligg^e  ability  of  inversion  to  estimate  a 


on  of  transmissivity.  The  analysis  requires  calc  jI 


,n  of  the  sensitivity  of  each  observation  to  the 


transmissivity  of  eachBell  in  the  model.  Incorporation  of  an  adjoint  states  based  sensitivity  calculation  into 
Dquri 


the  MOIlFLOW-2000qbrogram  significantly  speeds  thetletermination  of  the  sensitivity  matrix.  By  using 
tmnf^fJdKinQiilar  vhlue rlacnm^riKilinfi  we  hava  shnw3  for  our  hypothetical  example  that  multiple 
measurements  of  drawdown  at  a  single  location  have  similar  information  contents  with  respect  to 
transmissivity.  Resolution  analysis  can  pro\fiae  a  valuable  tool  for  both  assessing  and  designing  an 
inversion  procedure.  Distance  in  Meters 
Figure  4  Resolution  for  1  m 
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